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I. INTRODUCTION 

Observing strongly correlated atomic quantum gases in situ and real time is quite an achieve- 
ment [[Ill2l[3l|4l|51|6l|Vl[8l. By controlling the trapping geometry one can effectively adjust 
the "degree" of dimensionality, by feeding in more particles one can approach the thermodynamic 
limit, a central concept of our macroscopic world, and by controlling the coupling constant one can 
switch between universality classes of physical systems. Maybe, all of this was once envisioned by 
the great minds who have conceived the very few exactly solvable models of many-body physics 
m [TOl [m [l2l [El [Ml, but witnessing the merger of expectation and experiment proves to be an 
exciting period, today. 

In the current article we have explored a particular aspect of this rich topic, by focusing on the 
quantum properties of the ground state containing only a few bosonic particles, i. e. = 2, 3, 4 
inside a harmonic, one-dimensional trap. This situation is akin to the atomic or nuclear physics 
hmit of a condensed matter system. There the implied granularity of fermionic matter appears as a 
shell structure in the energy configuration or on energy surfaces through the appearance of magic 
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quantum numbers. This was unheard of in the field of uncharged gaseous matter until recently 
with the experimental use of neutral, repulsive, bosonic atoms [fT5l . 

A striking example here is given by the one-dimensional (ID) Bose gas. Its closed-form solu- 
tion has been given for the homogeneous system using Bethe's ansatz in the pioneering work of 
Lieb [11], which focused on the thermodynamic limit (iV ^ oo at fixed density). These results 
have been extended to include finite particle numbers [[T6l and effects due to a slowly varying 
trapping potential [fTTlfTSll . in the sense that the thermodynamic-hmit results for the homogeneous 
system still hold locally even in the presence of an inhomogeneity. However, for small atom num- 
bers A^, the trapped system can be solved in a numerically exact fashion, without resort to such 
a local-density approximation [fT9l l20l |2T1 |22l |23l. This is the starting point of our paper, which 
aims at a detailed comparison of the exact correlation functions of N trapped bosons with those 
obtained in a finite-size homogeneous system under a local-density approximation. That way we 
map out intrinsic confinement effects and discuss the validity of the local-density approximation 
for small atom numbers. 

Following this motivation, we will present the basic model of particles trapped in one di- 
mension in section [11} Next, we introduce our physical observables and computational methods 



in section III The homogeneous limit of this system is the Lieb-Liniger model, which represents 
our benchmark. Its basic notions will be reviewed briefly in section |IV] and used in section |V] to 
discuss our numerical results and analytical modeling of a few harmonically trapped particles. 



II. MODEL 



Let us consider a gas of one-dimensional bosons with repulsive, short range interactions in 
a harmonic trap Il20ll23ll24ll25]| . Then, the dynamics is given by the dimensionless Hamiltonian 

N N 

j=i j<i=i 
where we have measured energy in units of hu, length in multiples of the harmonic oscil- 
lator length aQ = y^h/muj and used the short hand notation dj = d/dxj. For exam- 
ple, such an effective one-dimensional description can be obtained by starting with real three- 
dimensional bosonic atoms of mass m in a strongly anisotropic external trapping potential 
V{r) = ^mtu'^x'^ + ^muj'j_(y'^ + z^). If the transverse level spacing is much larger than in the 
axial direction j3 = uj^/uj ^ 1, we can integrate out two dimensions by assuming that the two- 
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dimensional (yz) -subsystem only occupies the ground state. This procedure leads to the quasi 
one-dimensional couphng constant g = 2(3as/aQ, where denotes the s-wave scattering length 
of the bosons. 

Nowadays, this situation can be realized experimentally [15J and it is possible to investigate 
quantum correlations in situ. In particular, we are interested in the properties of the ground state 
of a few interacting bosons, i. e. = 2, 3, 4, . . . , and we will explore their quantum correlations. 
In a homogeneous system of length L = ao£ with a dimensionless scale i, the linear number 
density n = N/ £ is translation invariant and the qualitative behavior of the ground state strongly 
depends on the correlation parameter 'y = g/n. This was first described by Lieb and Liniger ffTTTl . 
In the thermodynamic limit limAr^_^oo N/i = n, it turns out that the state of the homogeneous gas 
of bosons is completely characterized by this parameter. It is customary to call bosons weakly cor- 
related for 7 ^ 1 (Gross-Pitaevskii regime) and strongly correlated for 7 ^ 1 (Tonks-Girardeau 
regime). 

The local density approximation extends this description to weakly inhomogeneous systems 
under the assumption that the variation of the ground state follows parametrically the spatial vari- 
ation of the single particle density. In the following, this assertion will be probed by explicitly 
constructing the A^-body wave function with the Bethe ansatz. We will analyze the behavior of the 
state over the whole range of interaction strengths for an increasing particle number {N = 2,3,4) 
in order to investigate the transition towards the thermodynamic hmit. The correlation functions 
of the ground state strongly depend on the particle number and this effect is most significant for 
few bosons. This analysis will provide a profound understanding of the inhomogeneous system, 
which can only be solved numerically otherwise. 



Assuming we have complete knowledge of the symmetrized and normalized A^-particle wave 
function '^{x) = '^{xi, . . . ,xn), then we need to extract relevant information about its behavior 
in terms of experimentally accessible observables [|26ll27l . Most relevant for this purpose are the 
number density n{x) = Np{x), which is proportional to the single particle density 



III. OBSERVABLES AND COMPUTATIONAL METHODS 




(2) 
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the first order correlation function measuring phase coherence 

(1). . J dX2---dXN'^*{x,X2,...,XN)'^{y,X2,...,XN) 

'[x,y) = ■^=^= (3) 

Vp(x)p(?/) 

and the diagonal of the second order correlation function measuring density fluctuations 

(2)/ X N -1 J dX3---dXN\'^{x,y,X3,...,XN)\^ 

9^'{x,y) = —— f \ r \ • (4) 

N p{x)p{y) 
We will focus on the behavior of the second order correlation function as it is a more sensitive 
probe for quantum statistical correlations in the system. Theoretically, much attention has already 
been directed towards second order correlation functions (IHl |25l |28l |29l [30l [HI [32l [3l, but the 
discussions were mostly concerned with their properties in the thermodynamic limit. In contrast, 
our interest is directed towards the detailed behavior of the second order correlation function for 
few boson systems, which differs from the thermodynamic limit. 

The next subsection is devoted to a brief introduction of the Multi-Configuration Time- 
Dependent Hartree (MCTDH) method which is used to calculate the A^-body ground state of the 
Hamiltonian in the presence of a trap in ([T]). 



A. Multi-Configuration Time-Dependent Hartree method 

The numerically exact MCTDH method [l34l |35]| is a quantum-dynamics tool which has been 
applied successfully to systems of few identical bosons [|2Tll22ll36ll37ll38l[39l . It solves the time- 
dependent A^-body Schrodinger equation {ihdt — H)^(x, t) = 0, as an initial- value problem by 
expanding the solution 

^i/{x,t) = J]aj(t)$j(a;,t), (5) 
Jec 

in terms of direct (or Hartree) product states ^j{x,t) = YliLi^jii^i^^) ^^'^ summing over all 
admissible configurations C = {J = (ji, . . . , jAr)|l < ji < s}. In turn, the still unknown, best 
single -particle functions < j < s} are represented in a fixed primitive basis imple- 

mented on a grid and s denotes the maximum number of required basis functions. The permutation 
symmetry of^{x, t) is ensured by the correct symmetrization of the expansion coefficients aj(t). 

Using the Dirac-Frenkel variational principle, one can derive equations of motion for both 
aj{t), and (f)j{x,t) [|35l . Integrating this system of differential-equations allows us to obtain the 
time evolution of the system via This has the advantage that the basis {$j(t)|J G C} is 



6 



variationally optimal at each time t. Thus it can be kept relatively small, rendering the procedure 
very efficient. 

Although designed for time-dependent simulations, it is also possible to apply this approach 
to stationary states. This is done via the so-called relaxation method [|40ll . The key idea is to 
propagate an initial wave function = 0) by the non-unitary, imaginary time propaga- 

tor U{t) = e^^'^. As r — i> oo, any excited state contribution is exponentially suppressed with 
^-{Em-Eo)T 2eft with the ground state. In practice, one relies on a more sophisticated 

scheme termed improved relaxation [41 J, which is much more robust especially for excitations. 
Here, the energy E = (\E'|if|\E') is minimized with respect to both the coefficients aj and the 
orbitals (pj{x). The effective eigenvalue problems thus obtained are then solved iteratively by first 
solving for aj with fixed orbitals and then optimizing (f)j{x) by propagating them in imaginary 
time over a short period. That cycle will then be repeated. 

Applying this procedure, we obtain the exact wave function for the ground state and can sub- 
sequently calculate the correlation functions according to (|3]) and (|4]). A profound understanding 
of the qualitative as well as quantitative details of this correlation function is the main topic of 
this contribution. For this purpose we compare the results of the MCTDH method to Lieb-Liniger 
theory for the homogeneous Bose gas which will be combined with a local density approxima- 
tion. However, before doing so we briefly review the main concepts of the Bethe ansatz and 
Lieb-Liniger theory, which can be solved exactly. 

B. Lieb-Liniger theory and Bethe ansatz for the homogeneous system 

As we have discussed, it is only possible to obtain the eigenstates and eigenvalues of the trapped 
Hamiltonian in ([!]) with a significant computational effort. However if one can disregard the 
external trapping potential and supply the system with periodic boundary conditions instead, we 
obtain the Hamiltonian 

N N 

Hll = Y.-¥1+11 35{x^ - xi) . (6) 
j=i j<i=i 

The corresponding eigenvalue problem has been solved analytically by Lieb and Liniger 
[fm [T6ll using the Bethe ansatz ^ and they derived the ground state properties - even for the 
thermodynamic limit. Before we discuss correlation functions and their functional behavior, we 
will briefly recall the quintessential steps. 
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To solve for the ground state and its energy, we have to consider the eigenvalue equation 

Hll'^^x) = E^{x), (7) 

where the spatial region h {x = (xi, . . . , xn)\0 < Xj < i} and the totally symmetric wave 
function satisfies periodic boundary conditions. The dimensionless system length i determines the 
linear number density n = N/i for a given particle number A^. It still occurs in our dimensionless 
formulation of the Lieb-Liniger Hamiltonian because we want to allow for a straightforward com- 
parison of the ground states of the Hamiltonians in ([T]) and Q for equal particle numbers N and 
equal interaction strengths g. Furthermore, having the length £ as a free parameter we can tune the 
number density n such that the correlation parameter in the homogeneous system is equivalent to 
the correlation parameter at a certain position in the trapped system. 

A subspace of the whole configuration space is the spatial simplex that contains only ascending 
coordinate A^-tuples, i.e. TZ = {x = (xi, . . . , xn) \ < xi < X2 < ■ ■ ■ < xn < i}- this region, 
^ together with its periodic boundary conditions are equivalent to 

N 

{J2-ld])^{x) = E^{x), (8) 
- d,)^{x) _ = g^{x) _ , (9) 

Xj^l—Xj Xj^l—Xj 

^{0,X2,...,xn) = ^(a;2,...,X7v,^), (10) 



da:'^(x,X2, ...,Xn) 



= d^'^{x2, ...,xn,x) 

x=0 



(11) 

x=e 



Due to the required symmetry of the wave function under particle exchange, knowledge of ^(aj) 
in the region TZ is equivalent to knowing \&(a3) in all other regions of the configuration space. In 



order to solve ([8HTT]), one uses the ansatz 

^{x) = ape"=-^ (12) 

where the summation extends over all A^! elements V of the permutation group Sj^. If we de- 
note the wave vector of the N particles by A; = {ki, . . . , k^), then the permuted vector is kp = 
. . . , fcp(Ar)). For convenience we also introduce the scalar product kj^x = J2f=i ^v(j)^j- 
Immediately, one obtains for the ground state energy of the -particle system 

1 ^ 

E = iNn-'eBiN,^), esiN,^) = _ ^(A:/)^ . (13) 
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The remaining task consists of determining the wave vector k such that all boundary conditions 
are fulfilled. After minor algebra which is outlined in [fTTI . one obtains the following equations 
for the components of the wave vectors of the ground state 

N 

{k^+, - k,)i = J2{e,+u - e,:) + 27r for j = 1, . . . , iV - 1 , (14) 



1=1 



9r^ = 2 arctan 



(15) 



Furthermore we note that the ground state solution possesses reflection symmetry, which means 
that for every positive component kj there exists a negative counterpart —kj. 

Before the ground state wave function can be calculated, we need to clarify how the factors ap 
are defined, such that the boundary conditions are fulfilled. This is done by first setting ai = 1. If 
V takes k into kp, then this is achieved by subsequent transpositions. For each transposition, the 
amplitude acquires a factor — e*^'^ , if ki and kj are transposed and ki is to the left of kj. The product 
of all these factors is a-p. Thus, we end up with the following wave functions for N = 2,3,4 
particles 

j^^iiHi+S-ii+krnx) _j_ gi(6»3i+e32+fc3i23:) _ gi(6'32+f3i+^'2i+fc32ia') (17) 

\E'(Xl X2 X3 X4) = e**'12343; _ gi(6'21+fc2134£c) _ g«(6'32+fcl324a') _ gi(6'43+fcl243a:) 
_|_gi(e43+6l21+fc2143:c) _|_ gi(6l3i+e21+fc2314a') _|_ gj(93i+6'32 + fc31243;) _|_ gi(6'42+S32+fcl342a') 
_^i{d32+9'il+d21+k-i2Ux) _ gi(6'41+6'31+6'21+fc234i:c) _ g«(6'42+f31+932+fc3i42£c) 
_j_gi(6'4l+932+6'3i+6»21+fc324ia') _|_ g«(6l42+6'43+fcl423;c) _ gi{6'43+6'42+6'32+fcl4323;) 
_gi(6l4i+e43+021+fc24133:) _|_ gi(6'43+6'41+931+S21+fc243ia') _j_ gi(e41+f42+S31+^'32 + fc34123:) 
_gi(S42+641+6'32+6'31+6'21+fc342l3:) _ gi(6'41+642+943 + fc4123a:) _j_ gi(6»4i+e43+6'42+6'32+fc4132:c) 
_^g*{e42+6'41+e'43+f2l+fc4213a;) _ gj(e42+e43+041+6'31+e21+fc423ia;) 

_gi('?43+f4l+6»42+031+e32 + fc4312£c) _^ gi(e43+042+6'41+032+6'31+f21+fc432l3:) (Jg) 

in the region 0<xi<...<a;Ar<£. 

In the thermodynamic limit it is possible to switch from a discrete distribution of k values to a 
continuous distribution and the ground state energy can then be written as [iQl fTTTl 

E = lNn'^e{-f) (19) 
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where 6(7) is given in terms of the solutions of the Lieb-Liniger equations 

„3 rl 



e(7) 



Me, 7) 



7 



A3(7) 
1 1 



de/^(e,7)r 



(20) 
(21) 



Using the Hellmann-Feynman theorem p2l. we can directly obtain the second order correlation 
function along the diagonal from the ground state energy 

d 



(2)/ 



^7 



eij(Ar,7), 



(22) 



which can be measured experimentally [[H [6l |43l |44l |451 |46l. In the thermodynamic limit, this 
simplifies to g^'^\x,x) = g^'^\0,0) = e'(7). 

LOr 




FIG. 1: Second order correlation function 5*^^^(0,0) vs. 7 for a homogeneous setup. Results from the 
thermodynamic limit (solid line) are compared to - from bottom to top - results for = 2, 3, 4, 10, 100 
bosons (dash-dotted lines). 



IV. CORRELATION FUNCTIONS OF THE HOMOGENEOUS SYSTEM 

To get a feeling for the second order correlation function and its dependence on the particle 
number as well as the correlation strength 7, we first consider the homogeneous system, where 
this is straight forward. The results are depicted in figure [T] where we compare ^''^^•'(O, 0) in the 
thermodynamic limit to results that are obtained by using finite values of the particle number A^. 
As can easily be seen, the thermodynamic limit is approached very quickly and for = 100 
bosons the exact results obtained with the Bethe ansatz are almost indistinguishable from the 
thermodynamic limit. 
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The decrease of the second order correlation function for large values of 7 is due to an inter- 
action induced antibunching also known as fermionization of bosons. In the limit 7 ^ 00, this 
behavior was first predicted by Girardeau [flOll and is due to a one-to-one correspondence between 
impenetrable bosons and spinless fermions in one-dimensional systems. This can be understood 
by looking at the limiting values of the wave vector k and the phases 6rs as a function of the in- 
teraction parameter 7. At large values of 7, the wave functions of ( 16 ), ( [T7| ) and ( 18 ) approach the 
form of Slater determinants. Thus, the wave-function of bosonic particles must vanish when 
two particles are getting close as if they were fermions. 




FIG. 2: Wave vector fc^/vr vs. 7 for = 2 (solid lines), A = 3 (dash-dotted lines), = 4 (dashed lines) 
and A = 10 (dotted lines). The functional form of the wave vectors is very similar for different particle 
numbers, there is only a different number of wave vector components for each particle number. 



In figure [2] we plot the results for the wave vectors for an increasing particle number N = 
2, 3, 4, 10. For odd particle numbers the components of the wave vectors approach even multiples 
of 27r, which is exactly the result that one would obtain for non-interacting, spinless fermions 
subject to periodic boundary conditions. However for even particle numbers, the components of 
the wave vectors are also separated by a constant spacing of 2n if we go to strong interactions, but 
they approach odd multiples of vr. 

This does not correspond to the behavior of non-interacting, spinless fermions subject to pe- 
riodic boundary conditions. This breakdown of the mapping between impenetrable bosons and 
spinless fermions was already pointed out by Girardeau and originates from the periodic bound- 
ary conditions. Considering hard wall boundary conditions, then there are no restrictions on the 
particle number A^. Furthermore it can be noticed that the dependence of the wave vectors on the 
correlation parameter 7 is very similar for the different particle numbers. Hence the two compo- 
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FIG. 3: Anti-diagonal of the second order correlation function g^'^\x, —x) vs. position x for the homo- 
geneous system. Results for A'^ = 2 are depicted in subplot a) for = 3 in subplot b) and for = 4 
in subplot c). The lines correspond to values of the correlation parameter of 7 = 0.1 (solid line), 7 = 1 
(dashed-dotted line), 7 = 10 (dashed line), 7 = 100 (thin dotted hne) and 7 = 00 (thick isolated dots). 

nents of the wave vectors with smallest magnitude for = 4 and A^ = 10 are indistinguishable 
from the solid line which depicts the components of the wave vector for A^ = 2. 

Apart from the behavior of the wave vectors themselves, a closer look on the wave function 
reveals that the factors ap appearing in (12) approach either +1 or —1 for 7 00. Hence the 
wave function approaches the limit of a Slater determinant as introduced by Girardeau. A simple 
way of understanding these behaviors is to consider the defining equations for the wave vectors 



(23) 



( 14 1 and ( 15 ) in the limit 7 00. Using a linear approximation for the arctan, one obtains 
- k,)i = 2rr-^ + 0(1/7^) , Oi, = ^^^^^ + 0(1/7^) 

which exactly yield the previously explained behavior for 7 — > cxo. 

In the remainder of this section we want to analyze the spatial dependence of the second order 
correlation function for the homogeneous problem. Due to the translational symmetry of the 
system, the correlation functions only depend on the relative distance \x — y\. Hence, it suffices to 
just study the anti-diagonal of the correlation functions for y = —x. In particular, we will examine 
the particle numbers A^ = 2,3,4, because they exhibit finite number effects which vanish in the 
thermodynamic limit [|23l . In figure |3] we depict the anti-diagonal of the second order correlation 
function g^'^\x, —x) for A^ = 2, 3, and 4 particles, respectively. For an increasing value of the 
correlation parameter 7, we again notice the strong antibunching at the origin, which reduces 
lim^^o 9'"'^^ (0,0,7) = 1 — 1 /N for ideal bosons and lim^_^oo 9^'^^ (x, —x, 7) = . 

Furthermore the off-diagonal develops an oscillatory behavior which becomes most pro- 
nounced for large values of the correlation parameter. The oscillatory behavior strongly depends 
on the particle number. For A^ strongly interacting particles in the system we expect A^ — 1 maxima 
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of g^'^^ {x, —x) in the interval < x < L/2. This can be understood by recalling the behavior of the 
wave vectors. In the limit 7 — > oo their components are all equally separated by 2n and therefore 
exactly A^— 1 different combinations of the components exist in the calculation of g^"^^ (x, —x) when 
we insert ( 12) in These combinations range from 2n to (A^ — l)27r and lead to trigonometric 
functions with — 1 different frequencies, thereby producing the behavior that was described 
above. Finally, it is also interesting to note that the anti-diagonal of the second order correlation 
function has a kink at a; = which is due to the delta interaction of the particles. However, in the 
limit 7 — > 00 the kink vanishes and g^'^\x, —x) approaches a smooth behavior at the origin. 

In general, it is possible to write down the analytic form of any correlation function because 
we have knowledge of the complete many-particle wave function. Using the particular form of 



the wave function for A^ = 2, 3, 4, as given in ( 16), ( 17) and ( 18 ), we have calculated the second 
order correlation function according to (|4]) by properly using the total symmetry of the respective 
wave functions. However, the lengthy analytical results are not enlightening and we refrain from 



writing them down in detail. Nevertheless, studying them in the limit 7 ^ 00, using (23 ), we find 
the following simple result 

,'^'(.. -X) = ^ (1 - E l^-^- ■ + OHM . (24) 

where we have also extrapolated the results for A^ = 2, 3, 4 to general particle numbers A^. This 
is consistent with the known limit for non-interacting fermions, as can be seen by evaluating the 
geometric series 



sm'^{27rxN/i) 
m sin2(27rx/£) 



g^^\x, -x) = 1- z: . (25) 



This result was already obtained by M. Girardeau [110]. However he first considered the limit 
7 ^ oo and used the fermionic wave function for the evaluation of the second order correlation 
function. 



V. CORRELATION FUNCTIONS OF THE INHOMOGENEOUS SYSTEM 

With this understanding of the second order correlation function and its dependence on the 
particle number A^ and correlation strength 7, we are now in the position to interpret the trapped 
results. For N = 2 particles, exact results for the ground state wave function are known B7ll48ll 
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and for = 3, 4 particles we used the MCTDH method. The corresponding results for the spatial 
correlations and a fixed interaction strength of (yf = 10 are shown in figure |4j 

To establish a common ground for the comparison of the homogeneous and the trapped system, 
we first of all restrict ourselves to the values of the second order correlation function in the center 
of the trap g^'^\x = 0, 0). Tuning the interaction from g = 0.001 50, we cover the whole range 
of the correlation parameter 7 from the weakly interacting GP regime to the TG regime of strong 
interactions. In figure [5] we compare the results in the trapped system for = 2,3,4 particles 
(diamonds, squares and circles) to the homogeneous results (solid, dashed and dashed-dotted line) 
for the same parameters of g and N. The length of the homogeneous system i is chosen such 
that the constant number density n = N/ £ equals the number density of the trapped system in the 
center of the trap. In general we can see a good agreement of the results with slight deviations at 
the crossover from weakly to strongly interacting bosons around 7 = 1. 




FIG. 4: Second order correlation function g^'^\x,y) for harmonically trapped particles versus the two- 
particle coordinate x,y for N = 2 in subplot a), for = 3 in subplot b) and for = 4 in subplot c). The 
interaction strength is g = 10. 

Thus we conclude that the correlation parameter 7 remains a valid parameter for the description 
of the inhomogeneous system as well. It is therefore possible to use a position dependent correla- 
tion parameter 7(x) = g/n(x), which is the central hypothesis of the local density approximation 
(LDA). 

In further considerations, we will apply this position dependent correlation parameter. Thus it 
is necessary to be acquainted with its spatial behavior. Hence, we plot the number density n{x) for 

= 2, 3, 4 particles for various strengths of the interaction in figure [6j With an increasing inter- 
action strength the density of the ground state changes from a Gaussian shape at weak interactions 
to a broadened distribution with Friedel oscillations [|49| for strong interactions. 

From the previous discussion we learn that the second order correlation function in the center 
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LO 




FIG. 5: Second order correlation function ^'^^^(O, 0) in the center of the trap as a function of the correlation 
parameter 7. The data corresponds to interaction strengths ranging from g = 0.001 to g = 50. The results 
for = 2, 3, 4 trapped particles (diamonds, squares, circles) are compared to a homogeneous system (solid, 
dashed, dashed-dotted hne). 




FIG. 6: Number density n{x) for A'' = 2 in subplot a), for = 3 in subplot b) and for A'^ = 4 in subplot c). 
The data corresponds to interaction strengths of g = 0.1, g = 1 and g = 10 (solid, dashed, dashed-dotted 
Une). 

ofthetrap, ^(2)(o 0), 

can well be understood by looking at the corresponding correlation function 
in the homogeneous system for the same parameters of g, N and 7. In the next step we want to 
investigate the behavior of the diagonal of the second order correlation function, g^"^^ (x,x). For this 
purpose we take the exact values of the diagonal of the second order correlation function for the 
trapped system (obtained with an MCTDH calculation) and plot them as a function of the position 
dependent correlation parameter 7(x). Comparing these values with the diagonal behavior of the 
second order correlation function in the homogeneous case, presented in figures [T] and [5j basically 
corresponds to a local density approximation. As the density decreases by moving out of the center 
of the trap, increasing values of 7 correspond to increasing values of x. 
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FIG. 7: Diagonal of the second order correlation function g^'^\x, x) as a function of the correlation param- 
eter 7(x). The data corresponds to interaction strengths of g = 0.001, g = 0.01, g = 0.1, g = I and 
g = 10. The results for = 2, 3, 4 for the trapped system (diamonds, squares, circles) are compared to the 
results for the homogeneous case (solid, dashed, dashed-dotted line) over the whole range of the correlation 
parameter in subplot a). A magnification of the behavior of the second order correlation function around 
7 = 15 is shown in subplot b). 



In figure [t] we plot the diagonal of the second order correlation function g^'^\x, x) for = 
2,3,4 at interaction strengths oi g = 0.001, g = 0.01, g = 0.1, g = 1 and g = 10. Generally 
speaking, the results for = 2,3,4 for the trapped system (diamonds, squares, circles) agree 
rather well with the results of the Bethe ansatz for the homogeneous system in the weakly inter- 
acting regime. However the trapped system deviates significantly from the homogeneous results 
for large correlation parameters. Starting with values of the correlation parameter around 7 ~ 1, 
we can see oscillations of the trapped results around the homogeneous curve. This means that the 
local density approximation gets less suited to describe the physical content of the trapped system. 

The breakdown of the local density approximation can well be understood if one considers 
the behavior of the density in the trap. For large interaction strengths the density develops the 
previously mentioned Friedel-type oscillations which can not be described within a local density 
approximation. This translates to the second order correlation function where similar oscillations 
occur. 

Whereas for = 3 the Friedel type oscillations lead to a peak of the density in the center of 
the trap, the opposite is true for = 2, 4, where we have a dip in the center. In terms of the second 
order correlation function this leads to oscillations around the homogeneous result that either start 
below the homogeneous curve (A^ = 3) or above it (A^ = 2, 4), as one moves out of the center 
of the trap. This can be seen in subplot b) of figure |7] which magnifies the behavior of the second 
order correlation function around 7 = 15. 
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FIG. 8: Second order correlation function g^^\x,y) in a local density approximation combined with the 
Bethe ansatz for = 2 in subplot a), for = 3 in subplot b) and for = 4 in subplot c). The interaction 
strength is in each case g = 10. 

Apart from an oscillatory behavior for large interaction strengths, we have seen that the diago- 
nal of the trapped system can be understood if we combine the Bethe ansatz with a local density 
approximation, in the sense described above. Finally, we want to investigate to which extent the 
full behavior of g^'^\x, y) can be analyzed by the same means. 

In this last step of the comparison we calculate the value of g^'^\x,y) in the local density 
approximation in the following way. In analogy to the previous discussion, we first of all recall 
that the comparison of the homogeneous and the inhomogeneous system is made for the same 
values of the particle number N and the coupling constant g. In the original homogeneous system 
we have the translational symmetry and the properties of the correlation functions only depend 
on the relative distance |x — y|. In an inhomogeneous system the relative distance is not the only 
relevant property that determines the behavior of the second order correlation function g^'^\x, y). 
Instead we have to incorporate the spatial dependence of the density to arrive at a local density 
approximation. For the diagonal part of the correlation function we have already seen that this 
combination of the Bethe ansatz and the local density approximation leads to a good agreement 
with the results for the inhomogeneous system. We can extend this procedure to non-diagonal 
coordinate pairs (x, y) by choosing the density at the center of mass coordinate [x + y)/2 for the 
local density approximation. As previously, this density is given by the exact MCTDH results and 
in the next step we again adjust the dimensionless length / for the homogeneous system, such that 
we obtain the same density for the fixed particle number N . Thereafter we solve the Bethe ansatz 
for these parameters and extract the anti-diagonal value of the second order correlation function 
for a relative distance of \x — y\. This result is eventually used to represent g^'^\x,y) for the 
combination of the Bethe ansatz and the local density approximation. This procedure is repeated 
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for every coordinate pair (x, y) which is used to plot the second order correlation function. 

In this fashion we obtain the counterparts to the exact trapped results in figure|4]and depict them 
in figure [8] The most striking difference that can be noticed by comparing the homogeneous to 
the trapped results is the large dip in the off-diagonal for = 2, 3 in the homogeneous case with 
local density approximation. This is merely due to the fact that the periodic boundary conditions 
in the Bethe ansatz prevent one from going to large distances. For a larger particle number and 
consequently a larger number density this difference begins to be negligible in the region of inter- 
est, as can be seen in the plot for = 4. Apart from this artifact the general features are similar 
and we conclude that we can also understand the overall behavior of the second order correlation 
function in terms of a local density approximation and the Bethe ansatz. 

Having a closer look at the anti-diagonal in the center of the trap in figure |9] we get a clearer 
illustration of these conclusions. While the periodicity prevents one from matching the exact 
physical behavior for any x, at least the initial slope and the approximate shape for x < 1 are 
modeled well. 




FIG. 9: Anti-diagonal of the second order correlation function g^'^\x, —x) across the center of the trap for 
= 2 in subplot a), for = 3 in subplot b) and for = 4 in subplot c). The interaction strength is (7 = 10 
and the MCTDH results (solid line) are compared to the combination of a local density approximation and 
the Bethe ansatz (dashed-dotted line). 

VI. CONCLUSION 

We have examined the ground state correlations for repulsive, quasi one-dimensional bosons 
in a harmonic trap. In particular, we have focused on the few particle limit = 2,3,4, .. ., 
where exact numerical solutions of the many particle Schrodinger equation are available with the 
Multi-Configuration Hartree method. These numerical results for the inhomogeneous system are 
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modeled with the analytical solution of the homogeneous problem using the Bethe ansatz and the 
local density approximation. Tuning the interaction strength from the weakly correlated Gross- 
Pitaevskii- to the strongly correlated Tonks-Girardeau regime reveals finite number effects in the 
second order correlation function beyond the local density approximation. 
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